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ABSTRACT 

We explore the weak lensing effect by line-of-sight halos and sub-halos with a mass 
of M < 10'' Mq in QSO-galaxy strong lens systems with quadruple images in a con- 
cordant ACDM universe. Using a polynomially fitted non-linear power spectrum P{k) 
obtained from A^-body simulations that can resolve halos with a mass of M ~ 10^ M©, 
or structures with a comoving wavenumber of fc ~ 3 x 10^ /iMpc~^, we find that the 
ratio of magnification perturbation due to intervening halos to that of a primary lens 
is typically ~ 10 per cent and the predicted values agree well with the estimated 
values for 6 observed QSO-galaxy lens systems with quadruple images in the mid- 
infrared band without considering the effects of substructures inside a primary lens. 
We also find that the estimated amplitudes of convergence perturbation for the 6 
lenses increase with the source redshift as predicted by theoretical models. Using an 
extrapolated matter power spectrum, we demonstrate that small halos or sub-halos 
in the line-of-sight with a mass of M = 10'^ — IO^Mq, or structures with a comoving 
wavenumber of fc = 3 x 10^ — 10^ /iMpc~^ can significantly affect the magnification ra- 
tios of the lensed images. Flux ratio anomalies in QSO-galaxy strong lens systems offer 
us a unique probe into clustering property of mini-halos with a mass of M < 10^ Af©- 
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1 INTRODUCTION 

Gravitational lensing is one of the most powerful tools for 
directly probing the structure and the distribution of dark 
matter. The remarkable agreement between the predicted 
and the observed weak lensing effects by large-scale struc- 
tures or clusters provides independent and consistent esti- 
mates of clustering property of dark matter on cosmic scales 
> 10 /i~^Mpc. However, we do not fully understand the clus- 
tering property on scales below ~ 1 /i~^Mpc, which corre- 
spond to individual galaxy halos. Although the cold dark 
matter (CDM) model predicts a large population of mini- 
halos (< 10^ Mq), the observed number of dwarf galaxies 
in our galaxy seems too low in comparison with the pre- 
dicted value. The discrepancy may be alleviated by some 
baryonic process, such as suppression of star formation by 
background UV rad i ation in the r eioniz ation epoch (e.g., 
iBullock et ai] (|200(]| '). iBusha et al.l (|2010|)'). or tidal disrup- 
tion due to a galactic disk (|D'Onghia et al. Alterna- 



tively, the suppression of the number count might be asso- 
ciated with super WIMPs or warm dark matter which has a 
larger free-streaming length than CDM (|Hisano et aklfioO^ l. 

In order to probe the clustering property of dark mat- 
ter at mass scales of < 1 /i~^Mpc, strong QSO-galaxy lensing 
systems with quadruple images have been used in literature 
(Metcalf fc Madaull200ll : IChiball2002l ) . In fact, the flux ratios 
in some quadruply lensed QSOs disagree with the prediction 
of best-fit lens models with a potential whose fluctuation 
scale is larger than the separation between the lensed im- 
ages. Such a discrepancy called the "anomalous flux ratio" 
has been consid ered as an imprint of substructure inside 
a lensi ng galaxy ([Mao fc Schneider ll998l:lMetcalf fc Madaul 



2001; Metcalf ct al. 2004; Chiba et al. 2005; Suga i et al 
2007: McKean et al. 2007: l yiorc ct al,.^ 20OT; Minezaki et al 
20091 : iMacLeod et '51(20091 : IVegetti et al.ll2012l v 
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However, recent studies based on high resolution sim- 
ulations suggested that the predicted substructure popula- 
tion is too low to explain the observed anomalous flux ra- 
tios iMaccio fc Miranda|[200^: lAmara et al. I I2OO6I : IXu et"al] 
l2009l . I2OIOI : IChenI l2009l : IChen et al.ll201ll ). More detailed 
modeling of gravitational potential of the lens on scales com- 
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parable to or larger than the dis tance between the lensed im- 
ages might also improve the fit l|Wong et al However, 
the origin of the anomalous flux ratios in some quadruple 
image systems such as B 1422+231 and MG0414+0534 has 
been veiled in mystery (|Chiba et al.l l2005l : iMinezaki et alj 
|2009[ ). 

In addition to substructures in lensing galaxy, any 
intergalactic halos along the entire line-of-sight from the 
source to the observer can perturb the lensing potential. 
Therefo re, they ma y chan ge the flux ratios of the lensed 
images. IChen et al.1 (|2003| ) have found that the contribu- 
tion from intergalactic halos modeled as singular isother- 
mal spheres would be < 10% of that fr om substructures 
within the lensing halo. iMetcalj l|2005al ) performed ray- 
tracing simulations for intergalactic halos with a mass of 
10** Mq ^ M s; 10^ Mq . Assumi ng that the halos have 
Navarro, Frenk & White (NFW) (.Navarro et al.lll997^ pro- 
files and the n umber density is given by the Press-Schechter 
mass function (jPress fc Schechteilll974 ). he found that four 
radio lensed QSOs that shows a strong cusp-caustic viola- 
tion are consistent with the predicted values without any 
substructures in the lensing galaxy. Assuming that halo pro- 
files are modeled as singular isothermal spheres and the 
num ber density is given by the Sheth-Tormen mass f unc- 
tion (ISheth fc Tormenll2002i '). lMiranda fc Maccid (|20o3) ob- 
tained a similar conclusion for three radio and two opti- 
cal/IR lensed QSOs. Using a A''-body simulation that can 
resolve halos with a mass of > 10* h~^MQ, and halos with a 
mass (10® Mq ^ M ^ 10* Mq) w hose number den sity obeys 
the Sheth-Tormen mass function, IXu et al.l (|2012h obtained 
a result that violation of the cusp-caustic relation caused 
by line-of-sight halos are comparable to (even larger than) 
those caused by intrinsic substructures though it depends 
sensitively on the halo profile. 

In order to estimate the magnification perturbation due 
to intervening halos more precisely, it is important to take 
into account various effects that have been overlooked in lit- 
erature. Firstly, if the shifts in relative positions of images 
and lens due to line-of-sight halos are too large, fitting a 
model with a smooth potential to the observed data becomes 
difficult since such a change is a consequence of a local effect. 
Moreover, even if the individual perturbing halo is not so 
massive, clustering halos could produce larger image shifts. 
Therefore, we need to incorporate the effects of clustering as 
well as the shifts of position of images and lens. Accuracy in 
observed positions of lensed images and lens would give an 
upper limit on the mass scale of perturbing halos. Secondly, 
in some lens systems, violation of the cusp-caustic relation 
might be caused by relatively massive faint s atellite galaxies 
in the neighborhood of the lensing galaxv (|McKean et al.l 
I2OO7I : IShin fc Evan3l2008l : iMacLeod et al.ll2009l ). Therefore, 
application of the cusp-caustic relation to generic lensed 
QSO systems may not be appropriate. Instead, we need 
to use other statistics to fit the model. Thirdly, the effects 
of massive line-of-sight halos should be subtracted off since 
they can contribute to low-order components in magnifica- 
tion tensor such as a constant convergence and an exter- 
nal shear in the lens model. Otherwise, we would estimate 
anomalies in the flux ratios systematically large because of 
double counting. 

In this paper, we explore the weak lensing effect due 
to line-of-sight halos in QSO-galaxy lensing systems taking 



these three effects into account and study how it will affect 
the flux ratios of lensed QSOs with quadruple images. To 
take into account of halo clustering, we use A^-body simula- 
tions to calculate the non-linear power spectrum of matter 
fluctuations down to mass scales of ~ 10^ H'^Mq. For sim- 
plicity, however, we do not put baryons in our A-body sim- 
ulations. Then we estimate the magnification perturbation 
using the obtained non-linear power spectrum and study 
wheather observed lensed QSO systems with quadruple im- 
ages are consistent with our model prediction. In section 
2, we describe magnification perturbation due to line-of- 
sight halos. In section 3, we derive analytic formulae for 
the power spectrum of convergence due to line-of-sight halos 
constrained from perturbations in image shifts. In section 4, 
we describe our A-body simulations for obtaining the non- 
linear power spectrum. In section 5, image shifts and magni- 
fication perturbation are investigated using a semi-analytic 
method developed in section 3. In section 6, we present re- 
sults on the fiux ratio anomalies using 6 samples of QSO- 
galaxy lensing systems with quadruple images observed in 
the mid infrared (MIR) band. In section 7, we conclude and 
discuss some relevant issues. In what follows, we assume 
a cosmology with a matter density ~ 0.272, a baryon 
density Q,i, — 0.046, a cosmological constant SIa = 0.728, 
the Hubble constant Ho — 70,km/s/Mpc, the spectrum 
index Us = 0.97, and the root-mean-square (rms) ampli- 
tude of matter fluctuations at 8h~^ Mpc, erg = 0.81, which 
are obtained from the observed CMB (WMAP 7yr result, 
(|jarosik et al.llioill )'). the baryon acoustic o scillations (Per- 
cival et al. 2010), and Ho l|Riess et al.ll2009l V 

2 PERTURBATION OF MAGNIFICATION 

Suppose a QSO at redshift zs is lensed by a primary lens- 
ing galaxy at zl to produce multiple images Xi and less 
massive intergalactic halos (secondary lenses) perturb the 
QSO-galaxy lens system. In what follows, we assume that 
the size of a light source is sufficiently small in compari- 
son with the Einstein radius of the primary lens and those 
of perturbers in the line-of-sight. Choosing coordinates cen- 
tered at a primary lens, given the angular position of a point 
on the source 0y, the angular position of the source dx lensed 
by the primary lens and the intergalactic halos is approxi- 
mately given by the lens equation defined at multiple lens 
planes i = 1, 2, ■ • ■ , N (see Fig. [T]), 

JV 

Ds0y=Dsex-J2D^,SCX,{x^), (1) 

1=1 

where Ds and Di^s are the angular diameter distances be- 
tween an observer and the source, the ith lens plane and the 
source, respectively, and oci and Xi are the defiection angle 
caused by a perturber and the position vector at the ith lens 
plane, respectively, r^'s satisfy 

_/ Doi0., i = l 

where Doi and Dij are the angular diameter distances be- 
tween an observer and the first lens plane, and between the 
ith lens and the jth lens planes, respectively. We assume 
that the primary lens is placed at the Ith plane. In general, 
it is difficult to solve equations ([1} and ((2)| since the position 
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Figure 1. Schematic diagram of ray tracing tlirougli multiple lens 
planes. Unperturbed light ray (dashed lines) starts from a source 
S and is deflected by the primary lens galaxy G at the primary 
lens (the Zth) plane and ends at O, a position of an observer. The 
light ray is perturbed by intergalactic halos in the line-of-sight 
(full curves), ai is the deflection angle at the ith lens plane for 
i = 1,2, - ■ ■ ,N. e is the shift of image position or that of the 
center of primary lens. 



Xi depends on any other positions Xj, j 7^ i. However, if the 
spatial derivatives of deflection angles are sufficiently small 
so that 



\Aa,\ 

\ai\ 



da.i 
dxi 



« 1, 



(3) 



is satisfied, then the hght ray approximately follows an un- 
perturbed geodesic 



DoiO^ - Di,iO.i{DL0:,), i>l, 



(4) 



and each deflection angle (Xi depends only on the position 
vector Xi that is independent of the other position vectors. 
This greatly simplifies the lens equations (O and © since 
each cti becomes independent each other. Then the inverse 
of the magnification tensor is approximately given by 



M~ 



d0y 



doii(ri) 



dai{xi) 



(5) 



In terms of convergence and shear 7^1 and 7^2 due to the 
primary lens at the position of a point-like lensed image Xi, 
the contribution from the primary lens can be written as 



dai{xi 



i^i + 7ii 

7i2 



7i2 
— 7il 



(6) 



In a similar manner, in terms of perturbations of conver- 
gence 5k and shear 571, ^72, the contribution from clustering 
dark matter in the line-of-sight can be written as 



5r = D-s^Y. 



doLi{xi 



Sk + 571 
572 



572 
Sk — Sji 



(7) 



The approximated lens equation at each image position is 
then 



(8) 



In the following, we assume that perturbations of fiux of 
images due to shifts of positions are sufficiently smaller than 



those due to distortion of the images. Then the perturbed 
magnification matrix is given by 

{Hi + 5fj.i)~^ = (1 - Ki — 7ii — <5«; — 571) 

X (1 - + 7,1 - <5k + 571) - (7i2 + (572)^, 

(9) 

where — {1 — Ki)^ — 7^^ — 7^2- A magnification contrast 
for image Xi is defined by = Sfii/fii. 

Up to linear order in fiiSit, HiSyi, and ^ii5"/2, the mag- 
nification contrast is approximated as 

2(1 - Ki)5K + 27^1^71 + 27i2(?72 /-.f^x 

which can be written as 

_ 2(1 - Ki)SK + 27,(571 

if magnification matrix for the primary lens is diagonalized 
(i.e., 7i2 = 0). 

We expect that line-of-sight structures that signifi- 
cantly perturb the fiuxes of images are relatively mas- 
sive halos, which add Sk > to the background conver- 
gence. Assuming matter fluctuations that are homogeneous 
and isotropic, the mean of the shear perturbation is van- 
ishing but the root-mean-squa re value is expected to be 
jBartelmann fc SchneideIj|2001^ 



(12) 



In what follows we assume that the shear perturbations sat- 
isfy 



Sk/V2 < S-y-i < Sk/V2, j = 1, 2. 



(13) 



The sign of magnification contrast depends on the cur- 
vature of the arrival time surface where the arrival time 
is stationary. If the arrival time is locally minimum, i.e., 
[1 — Ki)^ ~ 'jf > and 1 — K > 0, the density contrast 
satisfies 



5f (minima) > 



(2 - V2)yiSK 
(1 - K 



■Yi 



(14) 



since 1 — Ki > 7^ . As one can always choose local coordinates 
in which 7^ > 0, we have S^ > 0. If the arrival time is locally 
maximum, i.e., (1 — k^)^ — 7^ > and 1 — k < 0, we have 



(5f (maxima) < 



(2 - V2)^,5k 
' (1-K, 



It 



(15) 



since 1 — Ki < 7i, leading to 5f < 0. Thus strongly lensed 
images generated at a locally minimum /maximum point are 
magnified/ demagnified definitely by intervening massive ha- 
los. If a stationary point in the arrival time surface is a 
saddle one, i.e., (1 — Ki)^ — -yf < 0, we have 1 — < 7i and 
1 — Ki > — 7i as we assume 7^ > 0. Then the magnification 
contrast satisfies 



(1^^^ < ..(saddle) < .(^^^IMH. 



(i-K0^-7r 



(16) 



Therefore, the sign of magnification contrast cannot be de- 
termined definitely without additional conditions whereas 
the mean value is positive (S^) > for 1 — < and 
negative (5f ) < for 1 — Ki > 0. If the background con- 
vergence satisfies a condition 2|(1 — Ki)j > V^"/i, then the 
magnification contrast has a definite sign (5f < or 5f > 0. 
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Thus strongly lensed images generated at a saddle point tend 
to he demagnified/magmfied by intervening massive halos if 
1-Ki >0(<0). 

It is worthwhile to note that the mass-sheet degener- 
acy can be broken if intervening halos affect the fluxes of 
multiply lensed images significantly. For instance, suppose a 
transformation with a constant scalar A in the background 
convergence and shear 1 — Hi — > A(l — Ki) = f — and 
7i — ^ A7i = 7i, which preserves the positions of a lensed 
images of a point source by changing the position at the 
source plane as y — > Ay. Ifowever, the magnification con- 
trast 5f depends on A as 5f oc A"^ for |5f I < 1. Therefore, 
from observed 5^ , one would be able to put a constraint on 
A if &K and could be measured with sh ifts of positions 
of exte nded images due to intervening halos (jlnoue fc Chibal 
l2005al lbl). 

In order to quantify anomalies in flux ratios of a lens 
system with quadruple images, the cusp-caustic relation has 
been used in literature. In the positive cusp case, for close 
three adjacent bright images (A, B and C) with magnifica- 
tions piAjMs and ^c, and the opening angle 9 spanned by 
the center of two images in the ends and the lens center, the 
relation is 



written as 



Reus 



IMA + IJ.B + lJ.c\ 



0, 



(17) 



where \^j.a \ + ImsI + j^^cl oo and 6 — >■ 0. However, in prac- 
tice, none of observed quadruple image systems satisfy this 
asymptotic condition. For instance, the observed smallest 
opening angle is ~ 30°. Most radio or MIR quadruple 
image systems have even larger opening angles AS ~ 100°. 
Furthermore, the potential of the primary lens is sometimes 
not smooth. For example, luminous dwarf galaxies, or groups 
of galaxies in the neighborhood of quadruple images can sig- 
nificantly alter the flux ratios. In fact, some lens systems 
with anomalies in the fiux ratios may consist of multiple 
lenses. Although Rcusp is suitable for ideal systems with 
l/i^l + |/^s| + |/ic| CO and A9 ~ 0, it may not be suitable 
for most of observed quadruple-image systems. 

To circumvent this problem, we introduce a new esti- 
mator 



2(4 - Ns) 



E 



(5f -5'* (saddle, Kj < 1) 



(18) 



where 5^ denotes a magnification contrast for an image i 
with a positive parity or a negative parity with > 1 and 
5^ (saddle, Kj < 1) is a magnification contrast for an image j 
that has a negative parity with Kj < 1. Here Ns is the num- 
ber of saddles with k > 1. Magnification for an unperturbed 
system is given by a best-fit model based on positions of 
images and the center of the primary lens galaxy. Therefore, 
rj is not a directly observable quantity. Roughly speaking, rj 
corresponds to a mean magnification contrast per image due 
to clustering halos in the line-of-sight. Note that rj depends 
on only observed and modeled flux ratios provided that the 
magnification perturbations are sufficiently small. Here we 
put a negative sign before 5^ (saddle, < 1) because for 
systems with significant contribution from intervening ha- 
los, we expect demagnification for saddle points with k < 1, 
namely, {5^} < as we have seen. Suppose we have a set of 
images with two minima, A and C and one saddle B with 
Kb < 1- Then the estimator of flux-ratio anomalies can be 



(19) 



In terms of observed fluxes A, B, C and estimated unper- 
turbed fluxes A(i,B(i,Co, the estimator is approximately 
given by observed flux ratios, 

2 / y-iT-, \ 2n 



n 



[ ABo 
\AoB 



- 1 



+ (^-1 
\CqB 



(20) 



In a similar manner, for four-image system with two minima 
A and C and two saddles B and D with k < 1, the estimator 



7?2(A,B,C,D) 



+ 



8 
(SI 



(21) 



3 CONSTRAINED CONVERGENCE POWER 

We consider strong lens systems in which positions of mul- 
tiple images are well fit by the potential of the primary 
lens that consists of either a single or multiple lenses with a 
smooth potential though the image fiux ratios do not nec- 
essarily agree with the model prediction. In these systems, 
perturbation of image shifts due to other massive halos or 
voids in the neighborhoods of lensed images should be suf- 
ficiently small. This can be interpreted as an observational 
selection bias that no other dark massive halos or voids do 
not reside in the neighborhoods of line-of-sight of images 
since presence of these objects would otherwise perturb the 
positions of images and lenses significantly. Moreover, mod- 
eling the primary lens galaxies and neighboring groups or 
clusters would also induce an observational selection bias. 
The two-point correlation of matter density field in the line- 
of-sight constrained by these selection biases are determined 
by the best-fitting accuracy in position of images and lenses 
as follows. 

Suppose that positions of a pair of multiple images A 
and B with angular coordinates 6a and 9b separated by an 
angle 6ab are fit by a smooth lens model within an error e. 
This implies that total angular shifts S6 of image A and B 
due to intervening halos or voids in the line-of-sight should 
satisfy (see Fig.l) 



\se{eA)-5e{eA + 0AB)\ < 



(22) 



Assuming statistical isotropy and homogeneity for back- 
ground perturbations, equation p2|) gives 

(23) 



2{5r (0)) - 2{5e{o)se{eAB)) < e\ 



where {) denotes an ensemble average. 

Using Limber's approximation, a 2-point correlation 
function of astrometric shifts 58 for a pair of light rays sepa- 
rated by an angle 6 can be written in terms of matter power 
spectrum Ps( k;r), comoving d istance to the source rs, red- 
shift z{r) as l|Bartelmann fc Schneider., 2001i ) 

^se{e) = {5e{o)5e{e)) 



dk 



dr 



r — rs 
rs 



[1 + z(j 



— M/(fc; k,ut{r))P,{k- r)Mg{r)ke), 



(24) 
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where 



average of the estimator r;^(A,B,C,D) is given by 



air) 



rL{rs-r) 



r < vl 
r ^ TL 



(25) 



describes the trajectory of photons that pass through a pri- 
mary lens at comoving distance r = r_L, Jo is the zeroth- 
order Bessel function, and W (k; kcutir)) denotes the win- 
dow function in which modes with wavenumber k smaller 
than kcut{r) at comoving distance r — dr /2 < r < r + dr /2 
are significantly suppressed. We will discuss the property of 
W{k; kcut{r)) in detail in next section. From equations (I23|) . 
(|24p . and (|25p . one obtains the cutoff scale kcm as a func- 
tion of e and r. Because the accuracy in position fitting is 
generally far better than that of flux ratios in observations, 
it may still allow deviation in flux ratios due to constrained 
convergence and shear flelds in the line-of-sight. The con- 
strained 2-point correlation of convergence k as a function 
of a separation angle 6 is 



Ck(6I) = (5k(0)<5k(6»)) 



drr 



rs 



1 + z{r) 



f'°° dk 

X / —kW{k;k,^t{r;e))Ps{k;r)Jo{g{r)k9), 
Jo ^'^ 

(26) 

where g{r) is given by (|25|l . At small angular scales I 3> 
1, the constrained convergence power is given by a Hankel 
transform of equation p6|l 



P4i) = 2n j d9e^4e)Jo{ei). (27) 

In terms of obtained constrained convergence correlation 
function , one can estimate an ensemble average of the es- 
timator T]^ defined in the previous section, which measures 
anomaly in fiux ratios. For example, for three images with 
two minima A and C and one saddle B with ks < 1, using 
an approximation (lllf) . for \5^\ <C 1, an ensemble average of 
the estimator (I19II can be written as 



(Ja + Jb)c71{0) - 2JABi.eAB) 



+ (Js + Jc)ol{Q) - 2Jscd(fec) 



where 



and 



^?(4(1 



f + 27- ) 



(28) 



(29) 



(30) 



for i = A,B,C and 7i = (7^1 -|- 71^2)^^^- We have used a weU 
known fact that ^k(0) = 2^7„(S) and (571572) = (SkSjc,) ~ 
0, for a = 1,2 provided that background matter density 
fiuctuations are statistically homogeneous and isotropic. In 
a similar manner, for four-image system with two minima A 
and C and two saddles B and D with k < 1, an ensemble 



{ Ja + ./s)a2(0) - 2Jab^4Sab) 



+ {JC + JB)al{Q) - 2JcB£.n{ecB) 
+ (Ja + JD)<T2(0)-2JADe«(^AD) 

+ (Jc + JD)al{Q) - 2JccC«(6'cd) 



(31) 



where Ji and Jij are given by (|29[) and H30|) . Application of 
these statistics to mid-infrared lenses is presented in section 
6. 



4 NON-LINEAR POWER SPECTRUM 

In order to evaluate equations 1)24^ and (|26p . we need an 
accurate matter power spectrum at scales down to k~^ = 
0[1] kpc. However, analytical fitting formulae in liter- 
ature are not suit able for this purpo se. For example, the 
halo-fit model by ISmith et al.l (|2003l ) has been frequently 
used to evaluate the non- linear power spectrum P{k). How- 
ever, it was shown that this model underestimates the power 
spectrum by some tens percent than the latest cosmological 
simulation resul t s on small scales k > l/iMpc~^ (see e.g. 
iTakahashi et al.l \2Qli )'). Hence, in this study, we run cos- 
mological A'^-body simulations to investigate the non-linear 
power spectrum P{k) at galactic scales and we make a new 
fitting formula of P{k). For simplicity, however, we do not 
input baryon in our simulation. 

In our cosmological A''-body simulation, we use a cu- 
bic box with a comoving side length of 10/i~^Mpc with 



1024 and 512 coUisionless particles. We can check 



a numerical convergence of our simulation by comparing a 
low resolution simulation (A'^ = 512'^) with a high resolu- 
tion one (1024'^). The softening comoving length is fixed to 
be 2.5% of the mean particle separation, corresponding to 
0.25(0.5) ft"^ kpc for Nf, = 1024^(512^). The particle mass 
is 7.1 X 10*(5.7 X 10^)/i"^Mq for 1024^(512=^) coUisionless 
particles. Therefore, the minimum mass of halos resolved 
by our simulations is 1.4 x 10^ h~^MQ (corresponding to 20 
particles) . 

We use a si mulation cod e called Gadget2 

ijSpringel et al.l I2OOII : ISpringell |2005| ). We calculate the 
initial conditions of particles based on the second-order 
Lagr a ngian perturbation t heory (2LPT) (jCrocce et al.l 
I2OO6I : iNishimichi et all 120091') with the i nitial linear power 
spectrum obtained by lEisenstein fc Hul (|l999l ) . The initial 
redshift of our simulations is Zin = 99 and we dump the 
simulation results of the particle positions at 2: = — 4. 
We prepare two independent realizations for Np = 512'' 
at (2 = 0,0.35,0.7,1,1.5,2.2,3,4) and a single realiza- 
tion for A| = 1024^ at (2 = 0.35,0.7,1,1.5,2.2,3,4). 
In calculating the power spectrum P{k), we assign the 
particles on Ng — 1280'' grid using the cloud-in-cell (CIC) 
method to obtain density fluctuations. After performing 
the Fourier transform, we correct the window function 



of CIC as (5k 



2{Lk,/2Ng)] X 5k, where 5k 



is the density fluctuation in the Fourier space (Hockney 
& Eastwood 1988). In addition, to evaluate the power 
spectrum on small scales accurately, we fold the particle 
positions x into a smaller box by replacing x — >■ x%(L/2") 
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where the operation a%b gives the reminder of the division 
of a by 6 (e.g. Valageas & Nishimichi 2011). Then, the 
resolution becomes effectively 2" times finer. Here we use 
n = 2 and 4. Finally, we evaluate the power spectrum: 

where the summation is done from k—Ak/2 to fc+Afc/2 with 
a binwidth Afc, and A^^ is the number of mode in the bin. 
We use a logarithmic binwidth, A logj^g(fc//iMpc~^) = 0.05. 
For our polynomial fitting, we do not use the P{k) at small 
scales where the shot noise begin to dominate the signal. 
The Nyquist wave number determined by the mean parti- 
cle separation is fcNyq ~ {2n/L){Np/2), which corresponds 
to fcNyq = 320(160) /iMpc"^ for = 1024^(512^) parti- 
cles with L = 10/i~^Mpc. Hence, we can probe the density 
fluctuations on very small scales, k — 320 /iMpc~^. We have 
checked that the power spectra P{k) of our TV-body simula- 
tions agree with simulations with higher resolution in which 
we use finer simulation parameters of the time steps, force 
calculation, etc. within 2(6%) for k < 100( 320) /iMpc"^ 

W e also use simulation results by iTakahashi et al.l 
l|2012l ). They use the same codes as ours (Gadget2 and 2LTP 
initial condition) and employ 1024^ particles in the simu- 
lation boxes of L = 2000,800,320/1"^ Mpc and combined 
the P{k) on the different box sizes to cover a wide range 
of scales. They provide the P{k) up to fc = 30/iMpc~^ at 
z — — 10. Hence we use their result for k < 30/iMpc~^ 
and the present result for k > 30/iMpc~^. 

As one can clearly see in Fig. [5] the power spectra of 
our simulations at scales k > 30/iMpc~^ are significantly 
larger than the valu es predicted in the original halo-fit model 
l|Smith et al.ll2003l) . At scales k ~ 300/iMpc"^ which are 
relevant to the weak lensing effect of line-if-sight halos, our 
fitting formula, which improves the original halo-fit model 
down to fc ~ 300 hMpc~^ (see Appendix for details) predicts 
a factor of 2 — 3 enhancement in amplitude of the dimension- 
less power A(fc) = ^yA2{k) at rcdshifts < z < 2. As we 
shall see in the following sections, this enhancement plays 
an important role in explaining the origin of anomalies in 
fiux ratios. 

The suppression of power spectra P{k) predicted by 
our simulation at scales smaller than the Nyquist frequency 
is due to the lack of resolution. In fact, the power spec- 
tra P{k) in simulations with 1024^ particles are system- 
atically larger than those simulated with 512^ particles at 
k > 100 ZiMpc"^. In what follows, however, we use our new 
fitting formula to estimate the power spectra at very small 
scales k > 320 /iMpc~^. Although the accuracy is not guar- 
anteed, it may give a good approximation on very small 
scales since the power tends to increase as the resolution 
becomes higher. 

Assuming that the comoving size r of a density fluctu- 
ation is r ~ 7r/fc, the relation between a mass scale of clus- 
tering non-linear halos and the wavenumber can be given 
by 



M{k;z)^ —^/^ik}p{z), (33) 

where p{z) is a mean comoving matter density of the back- 
ground universe at redshift z. As shown in Fig. [S] for a 
given wavenumber k, the corresponding mass M increases 




k(h/Mpc) 

Figure 2. Plots of dimensionlcss power spectra, A-^(fc) = 
fc^P(fe)/(27r2), at z = 0, 0.35, 1 and 2.2 (from top to bottom). Our 
simulation results of L = 10/i~^Mpc box with 512^ particles are 
plotted as plus and cross symbols, corresponding to two indepen- 
dent realizations. Our simulation results of L = 10 fe~^M pc box 
with 1024'^ particles and those by ITakahashi et al. 1 1 |2012|) . which 
match our simulation results at fc = 30/iMpc~^ are plotted as 
filled circles and empty circles, respectively. Our fitting functions 

i see Appendix) an d the predicted values in the halo-fit model 
Smith et a"l]|2003l) are plotted as solid curves in red and gray, 
respectively. The vertical and dashed lines represent the Nyquist 
wave numbers AiNyq and the shot noises for = 10243 (left) and 
512^ (right), respectively. 



as the redshift z increases due to time evolution of the fluc- 
tuations. We find that a wavenumber that corresponds to 
the Nyquist frequency /cNyq = 320 hMpc~^ corresponds to 
M = 3x 10'' M0 at 2 = 0. fc = 1000 hMpc'^, 10000 /iMpc"^ 
corresponds to M = 1 x 10® Mq, 3 x 10^ Mq at z = 0, re- 
spectively. 



5 SEMI-ANALYTIC ESTIMATE 

Based on formalism developed in section 2 and 3, we can 
estimate the defiection and fiux change of strongly lensed 
images using the non-linear matter power spectrum obtained 
in section 4. Firstly, we calculate the rms of difference in the 
total angular shifts 50 between two images separated by 6, 

A5e{e) = [{{5e{e) - 5e{o) f)] (34) 

In order to study the dependence on the scale of matter 
fiuctuation in the line-of-sight, for simplicity, we adopt a 
"sharp k-space" window function, 

W^.(fc;fce.O = | ?' It'r^" (35) 
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ber k < k„ 



O[10^"^]/iMpc"^ is already taken into 




1 2 
log,o[^: Mpc/h] 

Figure 3. Approximated relation between mass M and a 
wavenumber k for non-linear matter fluctuations at redshift z = 
O(full), z = l(dashed), and z = 3(dot-dashed). 
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Figure 4. Rms value of difference in shifts of image position 
A59 as a function of a separation angle 9 for minimum wavenum- 
ber kmin = OhMpc"^ (dashed), lOOhMpc-i (full), SOO/iMpc"! 
(dot-dashed). Source and lens redshifts are 25 = 3, 2^ = 0.5, re- 
spectively. The upper limit of the wavenumber is assumed to be 
fcmax = 1000/iMpc-i. 



for which modes k < kmin = kcut are cut off in equation (|24p . 
As shown in Fig. [l] for kmm < 100 /iMpc~^, the difference 
in the image shift A50 becomes larger as the separation 
angle increases. However, for kmin ~ 500/iMpc~^, AS6 has 
a peak around 9 — 1.5 arcsec and it gradually decreases as 
9 increases due to a cutoff at small scales k > kmax- For a 
given separation angle, A58 gets smaller as kmin increases. 

Without any cutoff in the background fluctuations, it 
turns out that the difference in the shifts is AS9 ~ 0.06 
arcsec for separation angle = 1 arcsec, which is the typ- 
ical scale of the Einstein radius of a primary lens. This is 
significantly larger than the observed error e of relative posi- 
tions of light centers, which is of the order of 1 mas for most 
quadruple image lenses observed in the optical or radio band. 
This is because contribution from modes with wavenum- 



account in the fitted model, including neighboring clusters, 
groups, massive galaxy halos and luminous dwarf galaxy ha- 
los provided that the position of modeled images and the 
center of the primary lens agree with the observed values. 
In fact, as shown in Fig. [S] on angular scales smaller than 
1 arcsec, the largest contribution comes from modes on scales 
k ~ 10/iMpc~^ at which the one-halo term begins to dom- 
inate the matter power spectrum in the halo fitting model. 
For larger separation angles 8 > 1 arcsec, relative contribu- 
tion from smaller scale fluctuations decreases. If the residual 
difference in the position of lensed images is less than the 
rms value, then contribution from clustering halos on such 
scales should be negligible. As we can see in Fig. |6l if the 
cutoff scale kmin is constant as a function of redshift as we 
have assumed, then the dominant contribution comes from 
fluctuations near the lens plane at z ~ in the line-of- 
sight. This can be understood as a result of two effects :(i) 
The lensing weight function (r — rs){l + z{r)) in equation 
(19), is maximum around a half comoving distance to the 
source, and vanishes at boundaries at r = and r = rs- 
(ii) Because of convergence of photon trajectories toward 
a point-like source beyond the lens plane, the relevant co- 
moving scale of fluctuations becomes smaller, which leads to 
a suppression of the contribution. The latter effect (ii) be- 
comes more significant as a cutoff scale becomes larger and 
the lens redshift zl becomes smaller. 

Thus we interpret that the contribution from fiuctua- 
tions on scale fc ~ 10 hMpc~^ at ~ zl correspond to a pri- 
mary lens galaxy or luminous galaxies around the primary 
lens. If large-angle fluctuations of projected gravitational po- 
tential in the line-of-sight of multiply lensed images have 
been modeled as an external shear (m = 2) or low multipole 
terms (m = 3), which corresponds to contribution from clus- 
ters or groups of galaxies, we would only need to take into 
account small scale fluctuations with wavenumbers k > kcut 
where kcut is roughly equivalent to the scale of observable 
luminous dwarf galaxies at the lens plane. These small scale 
fluctuations may affect the flux ratios significantly if the 
total perturbations of convergence integrated in the line-of- 
sight is sufficiently large. 

In order to assess the effects of small scale fluctuations 
on the flux ratios, we consider rms of the self correlation 
of convergence o"re(0) = Vsk(0) as a function of source red- 
shift zs- We can see in Fig. [7] that contribution of small 



scale fluctuations ranging from fc„ 



0{m^)hMpc-' to 



kmax = 1000 /iMpc"^ would yields cr^(O) ~ 0.01 m conver- 
gence if the source redshift satisfles zs > 3. In other words, 
the surface density in small scale structures in line-of-sight 
is of the order of one percent of that of the primary lens. 
This may be too small to be of any importance. However, 
as we have seen in section 3, anomaly in the flux ratio is 
proportional to the magnification of the primary lens, i.e., 
rj ~ 2{jj,}aK.{0) where fi is the mean magnification of images 
provided that k ~ 7 ~ 0.5 and correlations between different 
images are negligible. For an image with modest magnifica- 
tion ~ 5, it would yield 10 per cent change in flux ratios, 
77 ~ 0.1 if zs ~ 3. Such a cha nge is sufficient to explain the 
order of observed anomalies (|Metcalj[2005al lbl). Moreover, 
if we take into account the correlation of convergence be- 
tween different images, anomaly in flux ratios can be more 
distinctive. As we can see in Fig. |8] for a separation angle 
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Figure 5. Contribution of each mode to A58 marginalized over 
redshift 2 as a function of wavenumber k. For separation angle 9 ~ 
1, major contribution to A59 comes from modes k ~ 10/iMpc~^. 
Plotted curves are for separation angle 6 = 0.5" (dashed), 1.0" 
(full), 1.5" (dot-dashed) with source redshift 2:5 = 3 and lens 
redshift = 0.5. We 10000 /iMpc-^ 



6 ~ 0.5 arcsec, the amplitude of 2-point correlation £,k{0) is 
still comparable to the self correlation cr^ (O). Therefore, we 
expect less significant anomaly in the flux ratios for systems 
with larger Einstein radius as long as accuracy in position 
fitting does not change. In other words, if such an anomaly 
is observed in the primary lens with large separation angles, 
the chance of significant perturbation to one of the lensed 
images is higher than the cases in which all the images are 
perturbed at the same time at similar levels. 

As is the case of shifts of image positions, the largest 
contribution to the amplitude of convergence (t„(0) comes 
from modes on scales k ~ lO^Mpc"^ at approximately a 
half distance to the source (Fig. [9|. If we consider a sharp 
k-space filter with kcut ~ kmin, then the small scale con- 
tribution to the convergence (t„ (0) is smaller than contribu- 
tion from modes with k ~ kmin ■ In other words, the largest 
contribution to the amplitude of convergence o-re(O) comes 
from fluctuations with wavenumber k ~ kmin at approxi- 
mately a half distance to the source. Therefore, we consider 
that influence of redshift dependence of filtering functions 
W{k;kcut{z)) on cr„(0) is small. 

So far we have considered the "sharp k-space" filter 
for cutting off the contribution from large scale fluctuations 
with k < kmin- However, in real setting, large scale fluctua- 
tions in the line-of-sight are removed in real space as com- 
ponents that yield an external shear and a constant conver- 
gence or distort the relative positions of the lensed images 
and the lens centroid. Therefore, contributions from modes 
k > kmin may not be negligible. In other words, more mas- 
sive objects placed near the line-of-sight could affect the flux 
ratios. In order to assess this effect, we consider a "Gaus- 
sian" fllter that is defined as an integration of the Gaussian 
distribution function as 



Wg{k;k, 



1 + Erf 



logi(,(fc/fcei,t) 



(36) 



V^/21oglo(l + <?). 
where Erf is the error function and q describes the width 
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Figure 6. Contribution of fluctuations on different planes at red- 
shift 2 to A59 marginalized over wavenumber k. Plotted curves 
are for minimum wavenumber kmin — 

20hMpc-l (dashed), 
lOOhMpc"^ (full), and 500/iMpc-l (dot-dashed). Thick and thin 
curves correspond to lens redshifts 2^ = 0.5 and zj^ = 2.99, re- 
spectively. We assume source redshift 25 = 3, separation angle 
9 = 1 arcsec, and the upper limit of the wavenumber kmax = 
1000 hMpc-i. 
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Figure 7. Self correlation of convergence as a function of source 
redshift 25. Plotted curves are for minimum wavenumber kmin = 
20/iMpc-l (dashed), lOO/iMpc"! (full), and 500/iMpc-i (dot- 
dashed). We assume lens redshift is z^ = 0.5 and the upper limit 
of the wavenumber is rCrnax 

= 1000 hMpc-i. 



of the filter A log^Q k ^ q. There should be an upper limit 
for q as the perturber in the line-of-sight is too massive, it 
becomes observable near the primary lens. As a reasonable 
guess, we consider two types of the "Gaussian" filter, q = 0.4 
and q = 0.9. Approximately 4 times massive objects are in- 
cluded for q = 0.4 and 20 times massive for q = 0.9 (Fig. llOf) . 
We choose two types of "UV" cutoff, kmax ~ 1000 /iMpc"^ 
and kmax ~ 10000 /iMpc~^. For a given maximally allowed 
shift e = AS6{6 = 1"), kcut is obtained from equations (24) 
and (25) for each filter. Then we calculate the 2-point cor- 
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Figure 8. 2-point correlation of convergence as a function of 
separation angle 8 between two images. The three curves are for 
k^i„ = 20 /iMpc-l (dashed), lOOhMpc-i (full), and 500 /iMpc^^ 
(dot-dashed). We assume lens redshift = 0.5, source red- 
shift zg = 3, and the upper limit of the wavcnumber kmax = 
1000 hMpc-^ 




Figure 9. Contour plots of d"^ (7^(0) / dzd Ink. We assume lens 
redshift zi^ = 0.5, source redshift 25 = 3, and the upper limit of 
the wavenumber kmax = 10000 /iMpc"'^. 



relation function of convergence ^k- We adopt separation 
angles between an image and a lens 6 — 0",0.5", and 1" 
as typical examples. As shown in Fig. Illl the differences 
between different types of filters are astonishingly small for 
6 = 0.5". For 9 = 0, 1", the relative difference in is at most 
35 per cent. Therefore, the effect of large scale fluctuations 
with k < kcut to rj is less than ~ 20 per cent. It should be 
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Figure 10. The "Gaussian" filters with q = 0.4 (dashed curve) 
q = 0.9 (dot-dashed curve) and the sharp k-space filter (full 
curve). We assume k^ut = lOOO/iMpc"^. 



noted, however, that for 6 = 0", amplitude of is system- 
atically reduced if contribution from large scale fluctuations 
is taken into account. This is because the contribution to 
the shift l\59 at 8 = 1" with respect to contribution from 
smaller fluctuations with k > kcut is larger than that to 
the 2-point correlation of convergence 5fc(0). If 5 > 0.9, we 
expect that the predicted rj will be much smaller. 



6 MIR QSO-GALAXY QUADRUPLE LENSES 

In what follows, we only use MIR data for flux of gravita- 
tionally lensed QSOs with quadruple images (see table 1). 
The number of lensed images used in analysis of flux ratios 
is denoted as A'^. We discard any data with small signal-to- 
noise ratio. Relative positions of macro-lensed images and 
the lensing galaxy are taken from the CfA-Arizona Space 
Telescope Lens Survey (CASTLES j3 

in the visible and near- 
IR bands except for H1413-I-117. The mean error in the rel- 
ative position between a lensed image and a position of the 
primary lens obtained from A*' samples is denoted as (e) . We 
assume that the errors of position of images or a lens are 
not correlated each other and obey Gaussian distributions. 
The mean separation angle between a lensed image and a 
position of a primary lens obtained from all the 4 images is 
denoted as (8) . We assume that the centroids of MIR images 
agree with those observed in the visible and near-IR bands. 

In this study, we do not use radio QSO lenses, which 
have been frequently used in literature because the finite 
source-size effect might be important in analyzing fiux-ratios 
perturbed by intervening clustering halos. The typical size of 
radio continuum emission region of QSO lenses is L ~ 10 pc. 
For lens systems with magnifications above jj, = 10, then the 
sizes of magnified images can be estimated as > 3 x 10 pc. 
Assuming that the size of the Einstein radius of is typically 
xe ~ 5kpc, a fractional difference in magnification with re- 
spect to a point-source is Sjj,/ fi ~ | ln{L/xE)\L/xE ~ 0.05 

^ See |http: / /www.cfa.harvard.edu/castles/| 
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Figure 11. 2-point correlation functions of convergence ^k{6) for 
the "Gaussian" filters with q = 0.4 (dashed curves), q = 0.9 (dot- 
dashed curve) and for the sharp k-space filter (full curve) as a 
function of the rms shift e in the relative position of a pair of 
images separated by 6. The upper and the lower panels show the 
plots for kmax = 1000 /iMpc"! and for k^ax = 10000 /iMpc"! , 
respectively. In each panel, the top, the middle, and the bottom 
three curves correspond to 6 = 0", 6 = 0.5", and 6 = 1", respec- 
tively. We assume zg = 3 and zj^ = 0.5. 



Einstein radius of the lensing galaxy is Sx/xe ~ 0(0.001) 
since the magnification perturbation due to shifts of Sx can 
be estimated as S^/^ ~ 5x/xe- 

As a fiducial model of these lenses, we adopt a sin- 
gular isothermal ellipsoid (SIE) plus an external shear 
(ES) (H ormann et al.|[l993 ). which can explain flat rotation 
curves. We use only relative positions of lensed quadruple 
images and the center of lensing galaxy for modeling. The 
parameters of the SIE plus ES model are the angular scale of 
the critical curve or the mass scale inside the critical curve 
6', the ellipticity e of the lens and its position angle 6^, the 
strength and the direction of the external shear {^,0^), the 
primary lens position on the lens plane {xq, J/o), and the im- 
age position {xi,yi). The angles 6e and 9-, are measured East 
of North expressed in the observer's coordinates. 

It should be emphasized that the observed MIR flux 
ratios are not used for making best-flt lens models. They 
are used for only estimating amplitudes of the expected flux 
anomalies rj. To find a set of best-fit parameters, we use a 
numerical code GRAVLENS developed by Keeton in or- 
der to implement the Xpos fitting of the positions, which 
have 10 degrees of freedom (8 for quadruple lensed images 
and 2 for the center of the primary lens). Because the SIE 
plus ES (SIE-ES) model has 9 degrees of freedom, resid- 
ual degree of freedom (dof) is 1. If Xpos/doi < 2 cannot be 
satisfied, we consider either a contribution from a luminous 
dwarf galaxy X modeled by an SIS with an Einstein radius 
bx in the neighborhood of the primary lens (SIE-ES-X) or 
introduce a large error for the position of the primary lens in 
order to satisfy a condition x^/dof < 2 (SIE-ES+). The lat- 
ter procedure may be verified in some lens systems because 
any unresolved luminous dwarf galaxies or inhomogeneous 
structures inside the galactic bulge of lens galaxy would shift 
the position of the center-of-light from the center of the lens 
potential. 

In what follows, we briefiy review our sample of QSO- 
galaxy lenses and the best-fit models. 



l|lnoue fc ChibalboOSal ) if a top-hat type source with an ap- 
parent size 1/ = 50 pc at the lens plane is placed at the center 
of an SIS. Even in more realistic cases, the order of the dif- 
ference would be the same as long as the potential has a form 
similar to an SIS. For instance, the correction term due to 
differential magnification is proportional to L/xe for an SIE 
lens because we have fi~^ ~ 1 — 0[xe/x] (|Kormann et al.l 
where x denotes the radial proper distance from the 
center of an SIE. Moreover, inclusion of substructure in the 
pr imary lens can boost th e perturbation by a factor of 2- 
3 l|Metcalf fc Amara|[20l3 ). Therefore, we expect ~ 10 per 
cent systematic change in the flux ratios. 

The size of continuum emission regions in the MIR band 
is typically much smaller than radio counterparts. In fact, 
the estimated source sizes in our sample of MIR lenses are 
L ~ 1 pc, which is significantly larger than the Einstein ra- 
dius of stars L ~ 0.01 pc. Therefore, our sample is free from 
the finite source-size effect and the microlensing effects due 
to stars. For 5 lenses in our sample, the point-like source 
approximation is valid at ~ 1 per cent level in flux ratios as 
long as magnification is not significantly large, i.e., fi < 10. 
Moreover, we note that the effect of differential magnifica- 
tion due to intervening halos is negligible as long as the or- 
der of shifts in the relative position of images divided by the 



6.1 B1422+231 

This is a cusp caustic lens that produces three colinear 
bright images A, B, and C with an image opening angle of 
77°. and a faint image D. The source is near a cusp in the 
astroid-shaped caustic. The obse rved MIR fiux ratios gives 
RcuBp = 0.20 (|Chiba et al.ll2005l ). This lens system is the 
first example that shows a vi olation of the cusp caustic rela- 
tion (|Mao fc Schneideill 19981 ). However, subsequent analysis 
revealed that the violation is not significant when marginal- 
ized over the opening an gle and the maximu m separation 
between the three images (iKeeton et al. l l2003l) . The redshift 



of the source zs ~ 3.62 ijKundic et al 



1993) is largest in 



our 6 samples and the primary lens is possibly an elliptic 
galaxy at zl = 0.34 (|Tonrvlll998l ). Although the posit ions of 
images and lens can be weU fit by the SIE-ES model l|Chibal 
|2002| ). the MIR fiux ratios between the image s A, B, and 
C are not consistent with the model prediction IChiba et al.l 
l|2005h . We have confirmed these results and found that in- 
clusion of m = 3 term with the external shear does not 
improve the fit to the MIR flux ratios. Therefore, it is likely 



See |http: / /redfive.rutgers.edu/~keeton/gravlens/ 1 
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that this lens system is perturbed by matter fluctuations 
on scales smaller than the separation between the images. 
Comparing the estimated flux ratios (table 2) to observed 
ones (table 1), it seems that image A, which is a maximum 
point in the time arrival surface is most likely to have been 
magnified by perturbers. 



6.2 MG0414+0534 

This is a fold caustic lens where a source is placed near 
an astroid-shaped caustic but not near a cusp in the caus- 
tic. The source at a redshift of zs = 2.63 9 is lensed by a 
foreg r ound elliptical galaxy at = 0.9584 (|Lawrence et al.l 
1 19951 : iTonrv fc Kochanekl ll999l V It consists of two close 
bright images Al and A2 separated by 0".415 and two 
fainter images B and C. We have found that the SIE-ES 
model does not give a good fit to the data. In order to im- 
prove the fit, we have considered a pos sible luminous satel- 
lite, object X jSchechter fc Mporelll993l) . as ha s been studied 
in previous lens models (e.g., Ros et all l|2000h '). The object 
X is modeled by an SIS at {xx, yx) = (0".857, 0".180) with 
an error of 0".01 as taken from CASTLES. The SIE-ES- 
X model yields a good fit to the positions of images and 
lens with = 0.003. However, the MIR flux ratio of Al 
to A2 is not consistent with the model prediction. The dis- 
crepancy remains ev en multipoles with m — 3 and m — 4 
terms are included (jMinezaki et al.ll2009l ). Therefore, it is 
likely that this lens system is perturbed by matter fluctua- 
tions on scales smaller than the separation between images 
Al and A2. Comparing the estimated flux ratios (table 2) 
to observed ones (table 1), it seems that image A2 which is 
a saddle point in the time arrival surface is most likely to 
have been demagnified by perturbers. 



6.3 H1413+117 

This is a "cross" lens in which quadruply lensed images have 
an approximate D4 (dihedral group with 4 rotational sym- 
metries) symmetry. It consists of 4 images, A, B, C, an d D. 
The source redshift is zs = 2.55 (|Magain et al.|[l98^ ') but 
the lens redshift is unknown. We use the d ata of the MIR 
flux r atios and the lens position observed by jMacLeod et al.l 
l2009l 'l. We have found that the SIE-ES model yields a good 
flt /d-of = 1.5 to the data of image and lens positions. 
However, this model gives a poor flt to the data of the 
flux ratios. The origin of anomalous flux ratios may be sub- 
structures inside the primary lens or cluster ing halos in the 
line-of -sight. To circumvent this problem, iMacLeod et al.l 
(|2009l ) added an SIS at the position of another galaxy G2 
lying at (-1. "87, 4". 14) with an error of 0.07" frorn imag e 
A, which corresponds to object #14 in lKneib et al.l l|l998l i. 
They have found that /dof ~ 1 for image and lens posi- 
tions, flux ratios, and weak pr iors for the lens parameters. 
ICoicoechea fc ShalvapinI (|2010l ). found that time delays be- 
tween images A-D are also consistent with this model SIE- 
ES-G2. Assuming a concordance cosmology with the Hub- 
ble constant Ho — 70kms^^Mpc~^ and density parameters 
Q.m = 0.3 and 0,^ ~ 0.7, they estimated the redshift of 
the primary lens as zi = 1.88t.o'^\- Although it is not clear 
whether G2 is the only component that would reproduce the 
observed flux rations, we adopt the SIE-ES-G2 model where 



the position of the center of an SIS is flxed to (— 1".87, 4". 14) 
with respect to image A. 



6.4 PG1115+080 

This is a fold caustic lens. The source at a redshift of 
Zs = 1.72 is l ensed by a foreground galaxy at zl ~ 0.31 
(iKundic et al ] [19971 It consists of two close bright images 
Al and A2 separated by 0".482, and two fainter images B 
and C. The data of position of images and lens are taken 
from CASTLES. We flnd that the SIE-ES model does not 
provide a good flt to the data unless the error in lens position 
is increased from 0.003" to 0.02" (SIE-ES+). Interestingly, 
the value is the same as the la error of lens position ob- 
tained by (|MacLeod et al.ll2009f) for H1413-I-117. This may 
be due to systematic problem s in determining th e position 
of the faint lens. As shown bv lChiba et al.l (|2005l ). the SIE- 
ES-I- model gives a good flt to the positions of the images 
and the lens and the MIR flux ratios. 



6.5 Q2237+0305 

This is the nearest lens in our sample with a "cross" conflg- 
uration of four images. The source is located at zs ~ 1.695 
an d the lens at zl = 0.03 94 (|Huchra et al.lll985l 'l. As shown 
bv lMinezakieTZI (Exi V the SIE-ES model gives a good 
flt to the positions of the images and the lens as well as the 
flux ratios. 



6.6 RXJ1131-1231 

This is a cusp caustic lens with a source at a reds hift of zs = 
1.72 lensed by a foreground galaxy at zl = 0.31 (|Sluse et al.l 
l2003h . Unfortunately, no data of fluxes due to the MIR 
continuum emission is available. Instead, we use the data 
of fluxes of [OIII] emission line from a narrow-l ine region 
(NLR) around the source QSO dSueai et al.ll2007l '). We flnd 
that the SIE-ES model does not provide a good flt to the 
data of positions of images and lens unless the error in the 
position of the primary lens is increased from 0.003" to 
0.017" (SIE-ES+). We flnd that the SIE-ES+ model gives a 
good flt to the positions of lensed images but the flt to the 
flux ratios is turned out to be not sufficiently good. Because 
the size of the NLR ~ 100 pc is significantly larger than the 
size of the MIR (near IR in rest frame) continuum emission 
region (~ 1 pc), we need a careful cons i derat ion on the finite 
source-size effect. In fact, ISugai et al.l (|2007l ') found a possi- 
ble imprint of an extended structure in the NLR region. The 
fractional contribution from the extended components can 
be ~ 20% for an aperture of 0".77 for lensed QSO images. 
The observed "bridge" between image A and image C and 
a shift of image B in the opposite direction to the critical 
curve suggest an asymmetric structure of the source. If this 
effect is taken into account, then the fit to the fiux ratios 
might be improved. This is because the differential magni- 
fication of images A and C can lead to a reduction in the 
fiux ratio I/^a/mcI as extended components of the images A 
and C are nearer to the caustic than the corresponding core 
components. 
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Table 1. Observed MIR Flux Ratios 



Lens 




ZS 


N 




Flux Ratio 




(") 


W (") 


Reference 


RXJ1131-1231(*) 


0.295 


0.658 


3 


A/B 


C/B 
1 ig+0.03 




0.017 


1.9 


1, 2 


Q2237+0305 


0.04 


1.695 


4 


B/A 
0.84 ± 0.05 


C/A 
0.46 ±0.02 


D/A 
0.87 ±0.05 


0.0046 


0.9 


1, 3 


PG1115+080 


0.31 


1.72 


2 


A2/A1 
0.93 ±0.06 






0.020 


1.2 


1, 4 


H1413+117 


1.88(**) 


2.55 


4 


B/A 
0.84 ±0.07 


C/A 
0.72 ±0.07 


D/A 
0.40 ± 0.06 


0.020 


0.6 


5 


MG0414+0534 


0.96 


2.639 


3 


A2/A1 
0.90 ± 0.04 


B/Al 
0.36 ±0.02 




0.0042 


1.2 


1, 3 


B1422+231 


0.34 


3.62 


3 


A/B 
0.94 ±0.05 


C/B 
0.57 ±0.06 




0.0042 


1.1 


1, 4 



Refercnccsil. CASTLES; 2. Sugai et al. 2007; 3. Minezaki et al. 2009; 4. Cliiba et al. 2005; 5. MacLeod 
et al. 2009 

Note: (*): [OIII] line flux ratios. (**): The lens redshift zi^ is obtained from a best-fit model using the 
observed positions of the images and the primary lens, the flux ratios, and the time-delays between the 
images assuming //q = 70km/s/Mpc. 



Table 2. Best-fit Model Parameters and Flux Ratios 



Model 


h' 
(") 


(a;o,yo) 
(") 


e 


fie 

(dog) 


7 


Or-.! 

(deg) 


(") 


dof 
(") 


Xpos 


Flux Ratio 






RXJ1131-1231 
SIE-ES+ 


1.83 


(2.039, 0.568) 


0.145 


-57.8 


0.120 


-81.8 




1 


1.3 


A/B 
1.66 


C/B 
0.909 




14.6 


Q2237+0305 
SIE-ES 


0.854 


(0.075, 0.939) 


0.371 


64.9 


0.015 


-46.8 




1 


0.004 


B/A 
0.887 


C/A 
0.447 


D/A 
0.825 


3.73 


PG1115+080 
SIE-ES+ 


1.14 


(-0.361, -1.342) 


0.156 


-83.0 


0.110 


51.8 




1 


1.0 


A2/A1 
0.912 






12.5 


H1413+117(*) 
SIE-ES-X 


0.561 


(-0.172, -0.561) 


0.204 


-14.5 


0.062 


55.7 


0.583 


2 


2.2 


B/A 
0.894 


C/A 
0.905 


D/A 
0.458 


5.24 


MG0414+0534(**) 
SIE-ES-X 


1.08 


(0.472,-1.277) 


0.232 


-82.1 


0.102 


53.8 


0.185 





0.003 


A2/A1 
1.039 


B/Al 
0.329 




13.1 


B1422+231 
SIE-ES 


0.755 


(-0.741,-0.658) 


0.309 


-56.6 


0.166 


-52.3 




1 


0.55 


A/B 
0.797 


C/B 
0.512 




4.91 



(*) Object X is modeled by an SIS whose position is fixed at {xx,yx) = (—1. "87, 4". 14). 

(**) Object X is modeled by an SIS whose position is fitted to {xx,yx) = (0. "857, 0". 180) with an error 

0".01. 

7 FLUX-RATIO ANOMALY In order to measure a possible contribution of clustering 

halos in the line-of-sight to the flux ratios, we calculate ry 
defined in equation (|18|) for our sample of 6 MIR lenses (5 
continuum and 1 line emission). We assume that the errors 
of the flux ratios and the positions of images and lenses obey 
the Gaussian statistics. We also assume that non-perturbed 
lens potentials are given by best-fit models using observed 
positions of images and lenses. We do not consider any con- 
tribution from substructures within a primary lens. Effects 
of image shifts on the fiux ratios due to intervening halos 
are assumed to be sufficiently small though this assumption 
might not be valid if the allowed shift is as large as e ~ 0.02". 
For simplicity, we use the sharp k-space filter for determin- 
ing the maximum scale of fiuctuations that can affect the 
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Figure 12. Plots of r; as a function of a source redshift zs for 
a sample of 6 MIR lenses (disk) and their predicted values for 
kmax = 320/iMpc~^ (circle), kmax = lOOOhMpc"^ (square), 
kmax = 10000 /iMpc~^ (solid triangle). Top: we assume e = (e), 
where (e) is the mean error in relative position of an image 
and a lens defined for A' images, (table 1). Middle: we assume 
e = 0.003", a typical value for a lensed image observed in the 
optical/IR band. Bottom: we assume a cut off fcjens due to lens 
modeling and e = 0.003". The error bars show the Icr errors in 
the observed MIR fluxes. 



flux ratios from errors e in the relative positions between an 
image and the center of the primary lens. To obtain the "IR" 
cutoff kmin , we use a mean separation angle {9} of an image 
and the center of the primary lens and a mean error (e) in 
relative positions of an image and the center of the primary 
lens obtained from quadruple images for each system. 

We consider three types of the "UV" cutoff, kmax = 
320, 1000, 10000 ^Mpc"\ kmax = 320/iMpc"^ corresponds 
to the Nyquist frequency /cNyq of our A'^— body simulation. 
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Figure 13. Plots of cutoff scales kcut and the source redshift 
zg for e in table (lower, blue) and for e = 0.003" (upper, violet) 
and fcjens (star, red). The "UV cut off" scales are assumed to 
be kmax = 320/iMpc~^ (circle), kmax = 1000 /iMpc""*^ (square), 
kmax = 10000 /iMpc"^ (solid triangle). 



0.015 




0.000 



Figure 14. Plots of approximated amplitude of convergence 
"■/{(O) ~ r]/2{fi) as a function of a source redshift zg for the 
observed MIR lenses (disk) and their predicted values for 
kmax = 320/iMpc~-'^ (circle), kmax = 1000 /iMpc"-"^ (square), 
kmax = 10000 /iMpc~^ (solid triangle) assuming kcut that 
corresponds to € = 0.003" and cut off ki^na due to lens modeling. 
The error bars show the Icr errors in the observed MIR fluxes. 
Assuming zj^ = 0.5, rigorous values of (Tk(0) are plotted as 
full curves for {kcut, kmax) = (300 hMpc"! , 320 /iMpc"!) 
(bottom), (750 /iMpc-\ 1000 /iMpc"^) (middle), and 
(1100 /iMpc-l, 10000 /iMpc-l) (top). 



For small scale fluctuations with wavenumber k > fcNyq, we 
extrapolate the power spectrum obtained in larger scales 
k < feNyq. It should be noted that the extrapolated power 
spectrum may be systematically larger/smaller than the 
correct value for k > fewyq. As shown in Fig. [T^] (top 
panel), we find that clustering halos with a mass scale of 
M < O[10^ Mq] or fc > 200/iMpc"^ in the line-of-sight are 
sufficient in explaining the observed anomalies in the fiux 
ratios. However, for Q2237-I-030 and PG1115+080, the pre- 
dicted anomalies seem too large unless the rms value of rj is 
at least comparative to the mean. The apparent discrepancy 
may be due to a possible systematic error in the position of 
the center of a faint lensing galaxy, which is significantly 
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larger than that of lensed images. A larger error e tends to 
give a larger eflFect on flux ratios. In order to see this sys- 
tematic effect, we also calculate rj assuming only errors in 
lensed image positions. Assuming that e = 0.003" for all the 
6 lenses ("constant e"), as shown in Fig. If 21 (middle panel), 
we find that the fit to the data of PGf f f5+080 is greatly 
improved. However, the fit to the data of Q2237-I-030 is not 
improved. This is because the lens redshift is exceptionally 
small {zl ~ 0.04) as compared to other 5 lenses in which 
zl ~ 0.3 — 2.0. As shown in Fig. If 31 the cut off scale for 
Q2237-I-230 is k^ut = 0[f- fO] /iMpc"\ This corresponds to 
mass scales of fO""" Mq, which are similar to mass scales 
of the primary lens. Therefore, the constraint from the shift 
of positions is not so stringent. However, we need to take 
into account the effect of lens modeling as well. Because a 
constant convergence and a constant shear of the primary 
lens are already taken into account in our model, we need 
to cut off modes that are equal to or larger than the size 
of the primary lens. For simplicity, we assume that modes 
with a half wavelength longer than twice the comoving ra- 
dius of the critical curve r{zL)6E of the primary lens are 
cut off, where r{zL) is the comoving distance to the lens at 
a redshift of zl- The corresponding cut off wavenumber is 
fciens = 2-K/Liens where Liens ~ Ar{zL)9E- Oe can be es- 
timated as a mean separation angle {0) between an image 
and the center of the primary for N images. As we can see in 
Fig. If 31 the cut off wavenumbers are fc;e„s ~ 3f 00 /iMpc~^ 
for Q2237-I-030 and fciens = O(fO^) /iMpc"^ for other 5 sys- 
tems. This means that the modeling effect is significant for 
systems in which the lens redshift is exceptionally small. 
Taking into account the modeling effect in addition to an 
assumption on the shifts of images and lens e = 0.003" po- 
sitions, we find that the expected rj for Q2237-I-030 is signif- 
icantly reduced (Fig. If 21 bottom panel). Furthermore, the 
predicted values for the other 5 lenses agree with the data 
at ~ 2 IT level without consideration of run-to-run variance 
of ?7. This result does not change even if we use conservative 
values e = (e) as shown in table f . Moreover, if we cut off 
the small scale modes k > ki^ns, the effect of differential 
magnification due to shifts of images is negligible in com- 
parison with the flux changes due to the weak lensing since 
the magnification perturbation due to the shifts of images is 
given by Sji/ ji ~ Sr/rE ~ 0(0. OOf), where the angular size 
of the lens is rE/rL ~ f" and the order of the image shift is 
Sr/rL ~ O(0.00f"). 

The result suggests that the flux ratio anomalies are 
caused by the weak lensing effect due to extragalactic halos 
with a mass M < lO'^ Mq in the line-of-sight. In order to see 
the source redshift zs dependence of the flux ratio anoma- 
lies, we plot the approximated rms amplitude of convergence 
due to intervening halos 



where (fi) is the mean magnification obtained from the ob- 
served images for a best-fit model (see table 2). As shown 
in Fig. If 41 the estimated (0) from the observed MIR flux 
ratios monotonically increases as the source redshift zs in- 
creases and it agrees well with theoretical prediction at ~ 2 cr 
level. It should be noted that the expected amplitude of con- 
vergence ffniO) for RXJff3f-f23f which shows a deviation 
at ~ 2 (J level might decrease if the finite source-size effect 



is taken int o acco unt as we discussed in section 6 (see also 
ISugai et al.l (|2007l )). 



8 CONCLUSIONS AND DISCUSSION 

We have studied the weak lensing effect by line-of-sight halos 
and sub-halos with a mass of M < f O'^ Mq in QSO-galaxy 
strong lens systems with quadruple images in a concor- 
dant ACDM universe. Using a polynomially fitted non-linear 
power spectrum P{k) obtained from A''-body simulations 
that can resolve halos with a mass of AI ~ fO^M©, or struc- 
tures with a comoving wavenumber k = 3.2 x fO^ /iMpc~^, 
we find that the ratio of magnification perturbation due 
to intervening halos to that of a primary lens is typically 
r; ~ O.f and the predicted values agree with the estimated 
values for 6 QSO-galaxy lens systems (continuum emission 
for 5 lenses, line emission from NLR for f lens) with quadru- 
ple images in the mid-infrared band without considering 
the effects of substructures inside the primary lens. The es- 
timated amplitudes of convergence perturbation for the 6 
lenses increase with the source redshift as predicted by our 
semi-analytical model. This feature strongly supports a hy- 
pothesis that the observed fiux ratio anomalies are caused by 
intervening halos rather than substructures associated with 
the primary lens. However, we do not exclude minor effects 
from substructures especially for systems with low lens red- 
shift Zl in which the weak lensing effect is small. Using an 
extrapolated matter power spectrum, we have demonstrated 
that small halos with a mass of M = f 0^ - f O^Mq can sig- 
nificantly affect the magnification ratios of lensed images. 

Instead of mass M, we have used comoving wavenumber 
k for parametrizing cut off scale of matter fluctuations due 
to intervening halos. We have considered two types of cut off, 
kcut and kiens ■ k^ut is determined from accuracy in positions 
of lensed images and the primary lens since intervening halos 
would induce shifts in relative positions of images, fcjens is 
given by the (effective) Einstein radius of the primary lens. 
Because large scale fluctuations are taken into account as a 
constant convergence and a constant shear in lens models, 
fluctuations that are larger than the Einstein radius should 
be neglected. Neglecting the shift of the center of a primary 
lens, we find that h^^s < kcut = O[f0^] /iMpc"^ for 5 MIR 
lenses and fccut ^ kiens = O[f0''] /iMpc~^ for f MIR lens in 
our sample. 

We have not used the cusp-caustic relation Rcusp in or- 
der to measure the strength of fiux ratio anomalies since 
most of our lens systems have either a complex structure 
(a luminous satellite) or a broad opening angle 9 > 30° . In- 
stead, we have devised a new statistic rj, to quantify the mag- 
nification perturbation. As we need a detailed lens model 
that fits the observed positions of images and lens, it may 
sounds less generic than using Rcusp- In fact, the mass- 
sheet degeneracy yields ambiguity in estimating the magni- 
fication perturbation. Different models with different radial 
profiles would certainly give different predictions. However, 
this is not a problem. As we have discussed, perturbations 
in convergence and shear 5k, S-y can be measured from ex- 
tended sources surrounding t he MIR continuum emitting 
region (|lnoue fc Chiball2005bD . From observed rj, 5k, and 
^7, we would be able to break the mass-sheet degeneracy. 
This means such an ambiguity can be removed by estimating 
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shifts of lensed images with spatial structures with respect 
to unperturbed ones. 

Because we have used a new statistic rj instead of Rcusp 
it is dif ficult to directly compare our result with previous 
studies l|Metcalj|2005al : IXu et al.ll2012l ) in which the effect 
of clustering halos is considered to be minor. However, as 
our numerically obtained non-linear power spectrum incor- 
porates all the effects of clustering halos and that of their 
substructures, our result indicates that clustering effect on 
mass scales of M < l(f Mq is much important than consid- 
ered in previous studies. 

In order to estimate the magnification perturbation con- 
strained from shifts of positions of images and lens, we have 
considered a "sharp k-space filter" for cutting off the fluc- 
tuations on large scales. If we use "Gaussian fllters" that 
are sufficiently smooth, variance in convergence can be sys- 
tematically decreased than using the "sharp k-space filter" . 
However, if we also consider a cut off due to modeling of a 
primary lens, such an effect may be negligible as large scale 
modes are taken into account as a constant convergence or 
shear. It should be noted that we have neglected effects of 
3 point or 4 point correlation of matter fluctuations, which 
may enhance the flux ratio anomalies. In order to incorpo- 
rate these effects and check validity of our approximation, 
we need to implement ray-tracing simulation based on A'^- 
body simulations. 

If we include effects of baryons, we naively expect fur- 
ther enhancement in magnification perturbation as baryon 
cooling would st eepen the gravitat i onal potential of halos 
at small scales (iRudd et all l2008l : ISemboloni et all I2OIII : 
Ivan Daalen et al.l l20Il''). Then our result would give a lower 
limit of the amplitude of perturbation in magnification ra- 
tios. However, feedback from supernovae or super massive 
black holes could suppre ss such a stee pening near the center 
of halo due to outfiows (|Booth fc Sch ave 2009). This might 
eventually suppress the magnification perturbation due to 
line-of-sight halos. Thus in order to improve our A'^-body 
simulations using only coUisionless dark matter particles, it 
is very important to incorporate baryonic physics down to 
mass scale of ~ 10'^ M0 or less. In other words, small scale 
baryonic physics which is relevant to galaxy formation might 
be gravitationally probed by the weak lensing effect in QSO- 
galaxy strong lensing system in the near IR or MIR band. 

Next generation telescopes such as the 
European Extreme Larg e Telescope (E-ELT) 
(jGilmozzi fc Spvromiliol l2007l'l or Thirty Meter Tele- 
scope (TUT) ijCrampton fc Ellerbroekl I2OO6I ) can be used 
to probe hundreds of such strong lens systems that are too 
faint for currently available largest telescopes to observe. 
They will provide us a unique probe into clustering property 



Survey with Subaru Telescope", by the Grand- in- Aid for 
the Global COE Program "Quest for Fundamental Prin- 
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MEXT Grant-in-Aid for Scientific Research on Innovative 
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surements of Images and Redshifts (SuMIRe)". Numerical 
computations were carried out on SRI6000 at YITP in Ky- 
oto University and Cray XT4 at Center for Computational 
Astrophysics, CfCA, of National Astronomical Observatory 
of Japan. 



APPENDIX A: FUNCTIONAL FORM OF 
FITTING FUNCTION 

In this Appendix, we provide the functional form of our 
fitting function of the non-linear matter power spectrum 
P{k). The fitting function can be used up to a wavenum- 
ber of = 320/iMpc~^, at ^ z ^ 4 for a con- 
cordant cosmological model with {^lm,^lt,QA,h,ns,crs) = 
(0.272, 0.046, 0.728, 0.70, 0.97, 0.81), which are obtaine d 
from the observed WMAP 7yr resuh ijjarosik et al.ll201lh . 
the baryon acoustic o scillations (Percival et al. 2010), and 
Ho (|Riess et al.ll2009l l. 

The non-linear power spectrum has been frequently cal- 
culated using the halo-fit model by Smith et al. (2003) (here- 
inafter, S03), which has 30 fitting parameters to fit the power 
spectrum obtained from their A'^-body simulations. However, 
as already pointed out by many authors (see. Introduction 
in Takahashi et al. 2012), the halo-fit model underestimates 
the power spectrum in comparison with the values from the 
latest simulations at small scales k > 0.1/iMpc~^. This is 
because the resolution of simulation in S03 is lower than that 
of the recent simulations. Recently, Takahashi et al. (2012) 
have provided an improved halo-fit model based on the 
original halo-fit model but re-calculated the fitting parame- 
ters to match their latest simulation results. However, their 
model can be used for only wavenumbers of < 30 /iMpc~^, 
which is not sufficient for our purpose. Our interest is in 
the galactic scale, which corresponds to fc > lOOftMpc^^. 
In this Appendix, we provide a fitting function that can 
be used to calculate the power spectrum at wavenumbers 
fc ^ 320 /iMpc~^. Our fitting function is based on the origi- 
nal halo-fit model (S03), but slightly change the parameters 
in their model to fit our simulation results. 

In the halo-fit model (S03), the dimensionless non- linear 
power spectrum, A^(fc) = k^P{k) / {2tv'^) consists of one- and 
two-halo terms: 



A\k)^A%{k) + Ai{k). 



(Al) 
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The first term is the two-halo term which dominates on large 
scales, while the second term is the one-halo term which 
dominates on small scales. We changed the one-halo term to 
fit our simulation results on small scales. In S03, there are 
four parameters a„,&„,c,i,7„ in the one-halo term. These 
four parameters are given by as polynomial functions of an 
effective spectrum index Uea and curvature C calculated by 
the input linear power spectrum (see Appendix of S03 for 
details). There are 17 coefficients in the polynomials, and 
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we will determine these 17 parameters to fit our simulation 
results. 

In order to obtain the fitting function, we use four sim- 
ulation results: three simulations with L = 10/iMpc~^ pre- 
sented in this paper (see section 4), and another from the 
work of Takahashi et al. (2012). For the four simulations, we 
adopt the same cosmological model as we have mentioned. 
The present three simulations with L = 10 ^Mpc~^ are used 
for fitting on small scales > 30 /iMpc^^ while the other one 
is used on large scales k < 30/iMpc~^. The fitting parame- 
ters are obtained by using the standard chi-squared analysis. 
The chi-square is defined as 



(A2) 



where Pmodei is the model prediction, and Pi,sim is the 
four simulation results labeled with an integer i = 
1 — 4. We sum up the powers at eight redshifts z = 
0, 0.35, 0.7, 1, 1.5, 2.2, 3, 43 and simply set the variance erf = 
^i^sim to give an equal weight for all the scales. 

For the three simulations with L — 10/iMpc~^, we sum 
up k to the Nyquist frequency /cNyq = 320(160) /iMpc~^ for 
Nl = 1024^(512^). On large scales k < 10/iMpc"\ the sam- 
ple variance among the three simulations is very large (over 
some tens percent) since our simulation box is very small 
(10/i~^Mpc on a side). Hence, we use the simulation results 
in the wavenumber where the sample variance is smaller 
than 16%, which corresponds to A; > 50^Mpc~^ for all the 
redshifts. 

We also use the simulation results from Takahashi et 
al. (2012) on large scales k < 30/iMpc~^. Using the same 
simulation codes as ours, they employed 1024^ particles in 
simulation boxes L = 2000, 800, 320 ^"^Mpc on a side and 
combined the P{k) from the different box sizes to cover a 
wide range of scales. They prepared 6(3) realizations for 
L = 320(800, 2000) h'^Mpc and gave the mean power spec- 
trum among the realizations up to fc = 30/iMpc~^. In or- 
der to reduce the shot noise effect, they did not use the 
simulation results at high wavenumbers k: the upper limits 
of the wavenumber are kmax ~ 30/iMpc~^ aX z = 0,0.35, 
fcma:r = 20/iMpc"^ at z = 0.7,1, fc,„ax = 10/iMpc"^ at 
z = 1.5,2.2 and kmax ~ 8/iMpc~^ at z = 3 in which the 
power spectrum is 10 times larger than the shot noise. 

Using the standard chi-squared analysis in Eo. (|A2[) . we 
find the best-fit parameters: 



logio an = 2.576 -f 2.263noff + 1.452neft + 0.6308nefi 

+ 0.1542ncft - 1.912C, 
logio &n = 2.062 + 1.034neff + 0.2651nfff - 3.677C, 



logi 



= 0.4449 + 1.743neff + 0.6772nfff + 0.06859C, 



7n 



= 0.2174 - 0.1366neff + 0.2418C. 



(A3) 



The other parameters such as On, /3n, /in, i^n, /(f^) are the 
same as in S03. The definitions of the effective spectrum 
index rieff and curvature C are given in the Appendix of 
S03. 

One can see in Fig[2]that our fitting formula agrees well 



^ We do not use an output at z = for one simulation of 
1024"^ in L = 10/t~^Mpc because of CPU time constraints. 



with our simulation results. In fact, the root-mean-square 
deviation of our best-fit model of A'^ (fc) from our simulation 
results is just 5.1% and the maximum deviation is 20% at 
fc ~ 50 /iMpc"^ due to the lack of available number of modes 
comparable to the simulation box size. 
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